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Three-dimensional Ising model in zero external field is exactly solved by operator algebras, similar 
to the Onsager's approach in two dimensions. The partition function of the simple cubic crystal 
imposed by the periodic boundary condition along both (010) and (001) directions and the screw 
boundary condition along the (100) direction is calculated rigorously. In the thermodynamic limit 
an integral replaces a sum in the formula of the partition function. When the z axis is chosen as the 
transfer matrix direction, a order-disorder transition in the infinite crystal occurs at a temperature 


T = T. determined by the condition: sinh A sinh 


EP = 1, where (Jı J2J) are the interaction 


energies in three directions, respectively. The analytical expressions for the internal energy and the 
specific heat are also given. It is also shown that the thermodynamic properties of 3D Ising model 
with Jı = J2 are connected to those in 2D Ising model with the interaction energies (J1J2p) by the 


relation fies E 


(es i 


PACS numbers: 05.50.+q, 64.60.-i, 75.10.-b 


I. INTRODUCTION 


The exact solution for three-dimensional (3D) Ising 
model has been one of the greatest challenges to the 
physics community for decades. In 1925, Ising presented 
the simple statistical model in order to study the order- 
disorder transition in ferromagnets [1]. Subsequently the 
so-called Ising model has been widely applied in con- 
densed matter physics. Unfortunately, one-dimensional 
Ising model has no phase transition at nonzero temper- 
ature. However, such systems could have a transition at 
nonzero temperature in higher dimensions [2]. In 1941, 
Kramers and Wannier located the critical point of two- 
dimensional (2D) Ising model at finite temperature by 
employing the dual transformation|3]. About two and 
a half years later Onsager solved exactly 2D Ising mod- 
el by using an algebraic approach [4] and calculated the 
thermodynamic properties. Contrary to the continuous 
internal energy, the specific heat becomes infinite at the 
transition temperature T = T, given by the condition: 
sinh mr sinh S = 1, where (J'J) are the interaction 
energies along two perpendicular directions in a plane, 
respectively. The partition function of 2D Ising model 
was also re-evaluated by a spinor analysis [5]. Up to now 
many 2D statistical systems have been exactly solved [6]. 

It is well-known that 3D Ising model has not been 
solved exactly yet due to its complexity. Because there 
is no dual transformation, the critical point of 3D Ising 
model cannot be fixed by such a symmetry. We also note 
that it is difficult to write the Hamiltonian along the 
third dimension of 3D Ising model with periodic bound- 
ary conditions in terms of the Onsager's operators. In 
addition, due to the existence of nonlocal rotation, 3D 
Ising model seems not to be also solved by the spinor 
analysis. Therefore, the key to solve 3D Ising model is 
to find the operator expression of the interaction along 
the third dimension. Here we introduce a set of opera- 


where z* = IIncothz = tanh "(e 


pm 


tors, which is similar to that in solving 2D Ising model 
[4]. Under suitable boundary conditions, 3D Ising mod- 
el with vanishing external field can be described by the 
operator algebras, and thus is solved exactly. 


II. THEORY 


Consider a simple cubic lattice with l layers, 
^ rows per layer, and m sites per row. Then 
the Hamiltonian of 3D Ising model is H = 
— an Cin aare + S207 Cj 41n + IRS R41) 
where o7;, = 1 is the spin on the site [ijk]. As- 
sume that rv; labels the spin configurations in the 
kth layer, we have 1 < vk < 2". As a result, 
the energy of a spin configuration of the crystal 
Ese = Ya Eri) + Dogar Balve) + Dar E(ie iei); 
where Ej(vj) and  E»(v) are the energies a- 
long two perpendicular directions in the kth 
layer, respectively, and E(vp,ve41) is the ener- 
gy between two adjacent layers. Now we define 


(ViV2) viv = 254 Qui (V2). = DVD (Vo) iwi = 
exp[-FA(v&)/(kgT) x  exp[-E»(vi)/(kpaT)] and 
(Ve)uw,, = expl[-E (vp, ve+1)/(kBT)). Here we 


use the periodic boundary conditions along both (010) 
and (001) directions and the screw boundary condi- 
tion along the (100) direction for simplicity [3] (see 
Fig. 1). So the spin configurations along X direc- 
tion in a layer can be described by the spin variables 
01,05,':: ,075,,. Because the probability of a spin 
configuration is proportional to exp[- E../(kpgT)] = 
(VMV2)u v, (V3) rave (ViVa) reve (V3) vous ub (V V2)uu (Vs)ua ; 
the partition function of 3D Ising model is 


` (ViV2)uv, (Va) 


(1) 


Z = Xava m Va Vv (VB ugs 
= tr(ViV2V3)!. 


FIG. 1: (Color online) The lattice structure in each layer of 
the simple cubic crystal. 


We note that Vi, V2 and V3 are 2""-dimensional matri- 
ces, and both V; and Və are diagonal. Following Ref. [4], 
we obtain 


V; = exp(Hi777, 070741) = exp(HiH;), 

V = exp( H2 Yi 0207.44.) = exp(H2 Hy), 2 

Vs = QsimhQH)" expl mor) O 
= [2sinh(2H)]"^/? exp(H* H,), 


where Hi = Ji/(kpgT), Ho = Ja/(kgT),H = J/(kpT), 
and H* = 4IncothH = tanh ! (e ?H). 

In order to diagonalize the transfer matrix V = 
Vi V2V3, following the Onsager’s famous work in two di- 
mensions, we first introduce the operators 


_ x ti. ee c z z 
Lia = Vt La,b = 049%a41%a4+2° °° C5 —10b (3) 


in spin space I along the X direction under the boundary 
conditions mentioned above. Here a,b = 1,2,--- ,2mn, 
07, oł and oz are the Pauli matrices at site a, respec- 


tively. Then we have L7 „= 1 and 
Lab mn — Latmn,b = -Q Lab = —LapQ (4) 


with Q = [[?", 02 = +1. It is obvious that the period 
of Lap is 2mn. We note that these operators La, are 
identical to P, in Ref. [4] except mn replaces n. 

H, and H, in the transfer matrix V can be expressed 
as 


He = Vigo leen an (5) 
H, = Dami Oz = Sol Loa: 


Following Onsager’s idea [4], we introduce the opera- 
tors 


1 2mn (a—b)rr 
ar = abc Lab COS ; 


^ 4mn B ( ur. 
mu 1 mn . (a—b)rm 
Br — —— 4mn aded Lab sin mn ? 
i 2mn " a—b)rr 
Ye = wig Daler (Laelt,2 — LaL) sin CHT 


(6) 


where z is an arbitrary index. Obviously, we have a_, = 


Qr, Bg = —Br, Bo = Bmn = 0, Yr = —^fyrs and Yo = 
Ymn = 0. Eqs. (6) can be rewritten as 


2mn 
Or = — Imr Dieci As COS mua: 
= 1 2mn sa rar 
Br — 2mn, » As SIN mnn’ (7) 
— 4 + rst 
Yr = 2mn Lal Gs SIN in? 
mn 
where A, = 2 lp Laats And G; = 


$3. ildbusbapaes — Jaabgara). According to 
the orthogonal properties of the coefficients, we obtain 


A, — Lesi [or cos Tet br sin Pez], (8) 
Gs = i$, ye sin t, 


From Eqs. (5)-(8), H, and H, have the expansions 


H, = A; = —2 SOS tus cos 77- — br sin 77-) 
— Q0 + Amn, (9) 


—1 
H, = — Áo = Qo + d Qr + Amn- 


Because Amn+s = —QA, = —A,Q and Gmn+s = 
—QG, = —G,Q, and combining with Eqs. (8), we have 


[1 + (-1)"Q]a. = [1 + (C1) QJB, = [1 + (-1)"Q]7r = 0. 
(10) 
When Q = 1, Q2r = Bor = Y?r = 0 while Q2p-] = 
Bor41 = Yar+1 = 0 if Q = —1. So we can investigate the 
algebra (8) with Q = 1 or -1 independently. However, we 
keep them together for convenience. In order to diago- 
nalize the transfer matrix V, we must first determine the 
commutation relations among the operators &r, 8, and 
Yr. Similar to those calculations in Ref. [4], we obtain 


[A Aj] = G5, G;] = 0, (11) 
[Gi Aj] = 2(Aj4i — A5-i)- 


Substituting Eqs. (8) into Eqs. (11), we arrive at 


lar, Br] = 2iyr, [Br Yr] = 2ia,, Mr, 6i] = 2ibr, (12) 


where r = 1,2,::- mn — 1, and all the other commu- 
tators vanish. Obviously, the algbra (12) is associated 
with the site r, and hence is local. Because a,, Br, and 
^j. obey the same commutation relations with —X,., —Y;., 
and —Z, in Ref. [4], we have the further relations 


ag = $(1— Q) = Ro, 

Onn = zll —(-1)""Q] = Rmn, 

r Br^fr = io, fr Or = ibr, (13) 
o2 2 y2 R? Rr, a = Rpa, = @&rRr, 

Br = RB, = Br Rr, Yr = F^. = oye Fo. 


We note that Am = paps = 

cd Ed Lia bg — Pond dani Lp (a-1)mpt(a- 1-5)m 
m 1 m n 

and Gsm = et Gps = phen 2 ail Lpt(a-1)ma x 
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FIG. 2: (Color online) Operator renormalization: schematic 
of £^ „in T, along the Y direction and L? , in T along the X 
direction. 


Ly (a—14s)m,z ie eee eee eee A where 
s =1,2,--- ,2n, and 
Áp, = Dit {aptam cos ore 
+Bp+(q—1)m sin [pig Dmjsr (14) 
Gps = DN Yp+(q—1)m Sin ipt en 


When m — p — 1, Eqs. (14) recover the results in t- 
wo dimensions [4]. It is obvious that A,; and Gp,j also 
satisfy the commutation relations (11). m p $p, 
[Apis Ap’, v= [Gp,5 Gy, y] = [Apis Gp, ve] — 

We have obtained the expressions of = and H, in 
terms of the operators a,, 8, and yr in the space T. In 
order to get the Hamiltonian in the third dimension, we 
project the operator algebra in the space T into the Y di- 
rection. Then we have m subspaces I';(p = 1,2,--- ,m), 
in which the operator algebra with period 2n is same with 
that in I’. In Tp, we define 


D —_ pT 

Lia m C 5--(a—1)m* 15 

EE. mg o? egt c? 9) 
a,b — *p4-(a—1)m* ptam pt+(b—2)m™ p+(b-1)m 


along the Y direction. Then we 
have Ap,s = ame eee and Gp,s = 
i Hm [Sota tont lade v Lr pt (a—1)m x 


te ae ee which also obey the same commutation 
relations (11) and (12), similar to Ap, and G,,,. Then 
the Hamiltonian Hy = » 5^ 4 Apa. 


Because [C ats: D5,,,] = 0 (see Fig. 2), we have 
[Ap,s, Ap,s] = 0, which leads to Ap, = A,,, due to their 
common local algebra (12). This is a renormalization 
of operator, which means that Aps and Ap s have same 
eigenfunctions and eigenvalues in I', or I' space. We note 
that V2 is the transfer matrix along Y direction, which 
must be calculated in I' rather than I’, space by mapping 
Ap,1 = Ap,1 in order to diagonalize total transfer matrix 


V. Therefore, we have 
Ay = psu 1 Apa = ae 1 Apa = = Ám 


-ao — 250 "(ap cos zz —B,.sin*t) (16) 
expe Omn- 


Here, we would like to mention that H; = —575 4 Apo = 
— Ao, which is same with that in (9). This means that 
when Jj = 0, the Hamiltonian of 2D Ising model is re- 
covered immediately. 

Because [Q, Hz] = [Q, Hj] = [Q. Hz] = [Q, V] = 0, V 
and Q can be simultaneously diagonalized on the same 
basis. In other words, the eigenvalue problem of V can 
be classified by the value +1 of Q. 

The transfer matrix V with Eqs. (9) and (16 


) becomes 


2 oHi Ai eHa Am o H* Ao 


V = [sinh(2H)| * e” 
= [2sinh(2 RES (H*—Hı-H2)ao (17) 
x lU elH* +Hı-(— 1)" H2lomn 


where 


U, = e 2H (ar cos {2 —6, sin Z£-) 
xe- 2H2 (ar cos 77 — f, sin TT) e2H* ar, 


In order to obtain the eigenvalues of the transfer ma- 
trix V, we first diagonalize U, by employing the general 
unitary transformation: 


ene ea (Or cos O,.-- B,. sin 0,.) U, 


x e^ Gr (Ar cos Or +B, sin Or) e- $NI — grow (18) 


Here 0, is an arbitrary constant and can be taken to be 
zero without loss of generality, and 


cosh £,. = D,., 
sinh £, cos Nr = Ar, 
sinh £, sin 7, = B, cosh(2a,.) 


tanh(2a,) = rs (19) 
— C, sinh(2a,.), 


where 


A, = cosh(2Hj) cosh(2H5) sinh(2H*) 

— sinh(2H1) cosh(2 H2) cosh(2H*) cos == 

— cosh(2H1) sinh(2H5) cosh(2H*) cos 17 
+ sinh(2H;) sinh(2H5) sinh(2H*) cos ŒV", 
B, = sinh(2Hi) cosh(2H5) cosh(2H*) sin 77- 
+ cosh(2H,) sinh(2H5) cosh(2H*) sin ™ 
+ sinh(2H;) sinh(2H2) sinh(2H*) sin CH7, 
C, = sinh(2Hi) cosh(2H35) sinh(2H*) sin 77- 
+ cosh(2H1) sinh(2H5) sinh(2H*) sin 77 

+ sinh(2H;) sinh(2H5) cosh(2H*) sin ŒH", 
D, = cosh(2Hi) cosh(2H5) cosh(2H*) 

— sinh(2H1) cosh(2 H2) sinh(2H*) cos mn 

— cosh(2 H1) sinh(2H5) sinh(2H*) cos *£ 

+ sinh (2H; ) sinh (2H) cosh(2H*) cos Ur. 


mn 


— 


— 


We note that D? + C2 — A? — B? = 1, which ensures 
that 3D Ising model can be solved exactly in the whole 


parameter space. When Hə = O(i.e.J2 = 0) and n = 1, 
or Hı = 0(i.e.J1 = 0) and m = 1, we have a, = H*. So 
Eqs. (19) recover the Onsager's results in 2D Ising model 
4]. 

Then the transfer matrix V has a diagonal form 


mn—1 


=Ï 34 —1 
eia  8üeYegl0Ta rArV e~ ITA arar 


mn—li : mn 
xe7 Urs1 ine = Dsinh(2H*)| 
xe UU -Hi-Ha)og t 32 Eran +[H* +H —(-1)™ Ho]amn 


(20) 


A. Transformation 1 


In order to explore the symmetries in 3D Ising model, 
we take the transformation 


ax = —Qr COS iin T + Br sin mno () 
Bt = ansin T + B, cos 2E, Y = -Yr 


It is easy to prove that až, 8; and 7; satisfy the same 
commutation relations with a,, 8, and yr. Then we have 


Hy = o$ 425m B cos P-D"T 4 ge sin (Ure) 
—(-1)" Gan: 

H, = —o$ — 29:2 [of cos “© — Bt sin ££] + o. 

(22) 

Obviously, such a transformation exchanges the interac- 
tion forms in (1,0, 0) and (0,0, 1) directions, but changes 
the interaction form in (0, 1,0) direction. Therefore, 3D 
Ising model has no a dual transformation, and the criti- 
cal point cannot be fixed by the Kramers and Wannier's 
approach [3]. 

'The transfer matrix can be expressed as 


V = [2 sinh(2H )] 7 eH141 eHz4m e-H* Ao 
= [Psinh(2H)]^£* e-- Ha -H")og (23) 
TI -U*e [Hi — (1) H2+H*]až mn, 
where 
Ut = glHiolg2Hiloz cos 72777 +87 sin CZ] 


Ed re * in DT 
xe-2H (or. cos £2 — 87 sin 45 


Following the procedure above, we can diagonalize the 
transfer matrix V, i.e. 


mn-—l ip mn—l * x mn—l * * 
eX 3-7 e22r—1 ara Ver 227-1 Q0 


xe- AC áma- 5 sinh(2H*)|* 
xe Hi Ha- H*)oeS YS 1 paž -[H4 —(—1)" Ho+H* Jax, 
(24) 
where 
sinh £, cos n = Až, tanh(2a;) = (25) 


— BE 
sinh €, sin y; = B; cosh(2a*) + C, sinh(2a*), 


and 
A* = sinh(2Hi)cosh(2H5) cosh(2H*) 
— cosh(2H,) cosh(2H5) sinh(2H*) cos 77- 
— sinh(2H,) sinh(2 H2) sinh(2H*) cos ** 
+ cosh(2H1) sinh(2H2) cosh(2H*) cos (marr 
B; = cosh(2H|i)cosh(2H;) sinh(2H*) sin = 
+ sinh(2H,) sinh(2H5) sinh(2H*) sin %7 
+ cosh(2H;) sinh(2H2) cosh(2H*) sin =U" 


We also have D2 + C2 — A”? — B? = 1. 


B. Transformation 2 


Let 
a, = —a,cos + 6, sin =, (26) 
B. = a, sin %7 ^. B, cos tT, yp = Irs 
we have 
Hy = 05 4272, [o]. cos BT — gr sin (6-7 
(1) amn; 
Hy = a0 T 2 xs Qj. F Pus 
H, = -ah — 290" al cos 77. — fj, sin ^£] 
(27) 
'The transfer matrix reads 
V = [2 sinh(2H)] 9e HiAi gHa Am g- H Ao 
= [2sinh(2H)] Y ei H2 - H*)o6 (28) 
x mE 1 U'el —(—1)” Hi-H3—(-1)" H*]ænn 
where 
Ul = e2H le. cos tenet re Bi sin (mtr ]o2Hao', 
A 


xe-2H' (an cos Z7 — B. sin m. 


Similarly, we have 


iE Tod gf mm-—l.,/,/ mm-—1l.,/,/! 
ex 2 Thr. ex 1 4 re Ve M 1 rar 
n—1 i 


xe- MT 2% = Dsinh(2H*)]^? 
x e Fi Ha- H*)og 3277 Era, -H[Ha- (71) (Hi tH" loin , 


(29) 


mn 


Here, 


tanh(2a’ — er. 
ni ( = WE (30) 
= B, cosh(2a;.) + C, sinh(2a;.), 


sinh £, cos n). = Al., 
sinh £, sin rj. 
and 
A'. = cosh(2H;) sinh(2H5) cosh(2H*) 
— cosh(2H1) cosh(2H2) sinh(2H*) cos 27 
— sinh(2H1) sinh(2H2) sinh(2H*) cos (mrs 
+ sinh(2H1) cosh(2H2) cosh(2H*) cos (n7 Dra , 
B, = cosh(2H|)cosh(2H5) sinh(2H*) sin == 
— sinh(2H1) cosh(2H2) cosh(2H*) sin (n—Drr 
+ sinh(2H;) sinh(2H2) sinh(2H*) sin C22", 


The identity D? + C2 — A? — B7? = 1 also holds. 


III. RESULTS 


Because o9,04,:::,O;4 have the common eigen- 
vectors xo with the corresponding eigenvalues 
Ao, /1,::* , Amn, from Eq. (20), we have Vw = Av, 
where 


Y = e^ D arar e7 Da” 3I ?r x, 
lnA— $mnIn[2sinh(2H)] + (H* — Hı — 
iyw lé, Ar + [Ay = (— 


H3)Ao 
1)? Ho + H*] Amn 

(31) 
At the critical point, we have £j = H* — Hı — Hə = 0 
[4]. This leads to a critical temperature T = Te given by 
the condition 


sinh(2H) sinh(2H, + 2H;) = 1. (32) 


If Hə = 0 or Hı = 0, we obtain the critical temperature in 
2D Ising model [3, 4]. We note that when H4 = H = H, 
the critical value He = J/(kpT.) = 0.30468893, which 
is larger than the conjectured value about 0.22 from the 
previous numerical simulations. 

We note that the thermodynamic properties of a large 
crystal are determined by the largest eigenvalue Amax of 
the transfer matrix V. Following Ref. [4], we have 


In Amax — mn In{2 sinh(2H)] 
& 6s toc E2z-1 for mn = 2L; 


33 

={& +3 +- + €on-1 (33) 
+ Hy — (-1)* H + H* for mn =2L +1. 

Here Ay = A3 — -- = Amn-1ı = 1, which are same 


with the eigenvalues of the operators X, in Ref. [4]. We 
note that these two results above can be combined due 
to €_, = €, and mn = 2[H, — (—1)™ H5 + H*]. So Eqs. 
(33) have the compact form 


In Amax — imn Inf? sinh(2H)] = 335554 £z 
zi pu) codi !Icosh(2H4) cosh(2H3) cosh(2H*) 
— sinh(2H;) cosh(2H») sinh(2H*) cos &7—D7 
— cosh(2H;) sinh(2H5) sinh(2H*) cos Grin 
+ sinh(2H;) sinh(2H3) cosh(2H*) cos @—=YCr— Vr), 
(34) 
In order to calculate the partition function per atom 
Noo = liMm noo naan) ron for the infinite crystal, we re- 
place the sum in Eq. (34) by the integral 


E )d, (35) 


—>00 


1 I 
In Ag; = 5 In[2 sinh(2H)] + 5; 4 lim 


where 


cosh Em (w) = D(w) = cosh(2Hi) cosh(2H2) cosh(2H*) 
— sinh(2H) cosh(2H2) sinh(2H*) cos w 
— cosh(2H;) sinh(2H2) sinh(2H*) cos(mw) 
+ sinh(2H;) sinh(2H2) cosh(2H*) cos|(m — 1)w]. 
(36) 


Similarly, the continuous A(w), A*(w), A’(w), B(w), 
B*(w), B'(w), C(w), &m(w), nw), n*(w), and m'(w) re- 
place the discrete Ar, A*, A’, Br, BX, BL, Cr, Er, Nr, ms 
and n}, respectively, by letting w = 77-. Here we empha- 
sis that when Hz = 0, or Hı = 0, Eq. (35) is nothing but 
the Onsager's famous result in the 2D case [4]. We also 
note that very different from the 2D case, the partition 
function of 3D Ising model is oscillatory with the system 
size m. Therefore, the conjectured values extrapolating 
to the infinite system seem to be unreliable. 

For a crystal of N — mnl, the free energy 


F-2U-TS--NkgTln As, (37) 


the internal energy 


U = F-T = NkgT? Ax 
On Ass O 1n Aes O 1n Ass 
—NkgT|Hi OH, + Ho OH» +H OH 1, 


(38) 


and the specific heat 


2 2 2 
C = 9p = NkplHi aia + He Gass + WG 
9? In Xoo 9? In Xoo 9? In As 
+2H, Ho 5H, die + 2H, H 8H.8H + 2H34H DHE | 
(39) 
Here, 
e —1 = limmo i cos 1" dw, 
oD 
oman = 4 limm—oo i, inh dw, 
oma — cosh(2H*) — + sinh(2H*) limmo fo cos ndw, 
2 
2 uM = 2 limm—oo jt sin? n* coth £,,dw, 
1 
O7 Indo _ : 
DH = + lim, 00 hu = Ieg (FE)? ] coth £j, du, 
2 
m. — up = 2sinh?(2H*)[- coth(2H*) lim, soo f cos ndw 
+ 2 limm+oo i sin 2 coth £z dio — 1], 
9? ids 


1 w (8A' | ƏD os 
OH, Os T z Mmo do "amt OH C087] coth Em), 


2 . 
mam ——il = sinh(2H*) 
aA* . 
. x Hittin 3 Jo necs — 2 cosh Êm cos n cos 7*), 
ô lns _ 1 * 
8H.BH = + sinh(2H ) 


dw ( OA oD 


H T 
X limino Jn sinh ém  \OH2 OH ; cos coth ém). 


We note that at the critical point, lim,,_,9 Em —^ 0. How- 
ever, lim,,-50 27- sinh Em — —cos7(0). Therefore, we 
can see from Eqs. (37) and (38) that at the critical point, 
the internal energy U is continuous while the specific heat 
C becomes infinite, similar to the 2D case. 

We consider the special case of J; = Jo, where the 
calculation of the thermodynamic functions can be sim- 


plified considerably. After integrating, Eq. (36) can be 
rewritten as 
cosh £s; (w) = cosh(2H1) cosh(2H5;;) (40) 


— sinh(2H,) sinh(2H3,,) cos w, 
where 


Hip = H* — Hy. (41) 


It is surprising that Eq. (40) is nothing but that in 
2D Ising model with the interaction energies (Ji, Jop) 
and Həp = 2%. Therefore, In As; — i In[2sinh(2H)J 
in three dimensions can be obtained from ln A2P — 
iIn[2sinh(2H5p)] in two dimensions by taking the trans- 
formation (41). In other words, the thermodynamic 
properties of 3D Ising model originate from those in 2D 
case. We can also see from Eq. (41) that both 2D and 3D 
Ising systems approach simultaneously the critical point, 
ie. Hžp = Hı and H* = 2Hj. It is expected that the 
scaling laws near the critical point in two dimensions also 
hold in three dimensions [6]. 


The energy U and the specific heat C of 2D Ising model 
with the quadratic symmetry (i.e. Hı = Həp) have been 
calculated analytically by Onsager and can be expressed 
in terms of the complete elliptic integrals [4]. The criti- 
cal exponent associated with the specific heat agp = 0. 
Because 3D Ising model with the simple cubic symmetry 
(i.e. Hı = Hə = H) can be mapped exactly into 2D one 
by Eq. (41), the expressions of U and C in three dimen- 
sions have similar forms with those in two dimensions. 
So the critical exponent a3p of the 3D Ising model is 
identical to agp, i.e. a3p = 0. According to the scaling 
laws dv = 2 — a and p +v = 2 — a [6], we have v3p = $ 
and uap = 4. 

Up to now, we have obtained the partition function 
per site and some physical quantities when the z axis is 
chosen as the transfer matrix direction. However, if the 
x(y) axis is parallel to the transfer matrix direction, the 
corresponding partition function per site can be achieved 
from Eqs. (35) and (36) by exchanging the interaction 
constants along the z(y) and z axes. Therefore, the total 
physical quantity in 3D Ising model, such as the free en- 
ergy, the internal energy, the specific heat, and etc., can 
be calculated by taking the average over three directions. 
We note that the average of a physical quantity naturally 
holds for 2D Ising model. 


IV. HIGH TEMPERATURE EXPANSIONS 


Now we calculate the high temperature expansions of 
the partition function per atom when J; = J2 = J. Ac- 
cording to the identity 


27 
1 In(2coshz — 2cosw’)dw’ = 2n, (42) 
0 


from Eqs. (35) and (36), we obtain 


InAs = sty fo Jo In(cosh?(2H) 
—sinh(2H)cosh(2H)|[cosw + cos(mw)] 

+sinh?(2H)cosh(2H)cos[(m — 1)w] 

—sinh(2H)cosw’ }dwdw’ 

= 3lncoshH + 3In(1 + k?) 


T PT 2k(1—k? 
tam Jo Jo In{1 REESE 


2 
qs cos[(m 1)w] 
. 2k(1—k?)? 


[cosw + cos(mw)] 


rrey cos’ j duodu' 
= 3lncoshH — 3k* — 62k6 — 298148 
—21024k19 — ... 


(43) 
where k — tanh H. Therefore, the partition function per 
atom in high temperatures is 


doo = 2cosh® H (1—3k*— 62k6 — 1036k? —20838k19 — . .. ). 
(44) 

We note that for periodic boundary conditions, the high 

temperature partition function per atom reads [7] 


AP. = 2 cosh? H (1-- 3k* + 22k° + 192k8 + 2046k19 4+ . . . ). 

(45) 
Obviously, the difference between A, and A2, comes from 
the screw boundary condition along the X direction (see 
Fig. 1). 


V. CONCLUSIONS 


We have exactly solved 3D Ising model by an algebraic 
approach. The critical temperature T;.(i = 1,2,3), at 
which an order transition occurs, is determined. At T;., 
the internal energy is continuous while the specific heat 
diverges. Obviously, the thermodynamic properties in 
three dimensions are highly correlated to those of 2D 
Ising system. When the interaction energy in the third 
dimension vanishes, the Onsager's exact solution of 2D 
Ising model is recovered immediately. This guarantees 
the correctness of the exact solution of 3D Ising model. 
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